/****************************************************************************************************************
*****************************************************************************************************************/
/*****************************AIM 3 INTIAL ANAYSIS***************************************************************/
/*******************STEP 1: EXAMING THE FERTILITY INTENTIONS VARIABLE*******************************************/
proc freq data=mydata.testaim2;
tables pregint*certain;
run;
proc freq data=mydata.LCAaim2;
tables pregint*certain;
tables intent*reactslf;
run;

/*proc freq data=mydata.LCAaim3;
tables intention intentcoll intention*intentcoll;
run;*/

/**********************STEP 1B: CREATE AND HARD ASSIGN LCA VARIABLE************************************************/
proc print data=postsx1; /***posterior probabilities from random sample*********/
run;

proc LCA data=mydata.LCAaim2 outpost=mydata.postprob1 outparam=par4 outstderr=err4 outseeds=seed4;
 Title 'LCA round 3 run menbmi categories';
	nclass 3;
	items sexprt2 sexrisk condom2 hltstat2 drink2 drug sti bmi3;
	categories 2 2 2 2 2 2 2 2;
	seed 123456789;
	id caseid;
	nstarts 50;
	maxiter 50000;
	rho prior= 1;
	weight weightvar;
	cluster panelvar;
run;


/********************STEP 1C: MERGING ORIGINAL DATASET WITH POSTERIOR PROBABILITY OUTPUT DATASET************************/

proc print data=mydata.postprob1 (obs=50); run; 
/*proc freq data=mydata.LCAaim2; tables kiddos; run;*/


proc sort data=mydata.LCAaim2;
by caseid;
proc sort data=mydata.postprob1;
by caseid;

data mydata.LCA1119aim3;
merge mydata.LCAaim2 (in=a) mydata.postprob1;
by caseid; if a;
run;

proc contents data=mydata.LCA1119aim3;run;

proc freq data=mydata.LCA1119aim3; tables best; run;
/****************************************CREATING PHP VARIABLE*****************************************************/

proc format;
value phpx 0="lowest risk"
		   1="substance users"
		   2="sexual risk takers";
value intentcollx 0="definite no"
	              1="definite yes"
				  2="Far or contingent intention"
				  3="indifferent or ambivalent";
label php = "preconception health"
      intentcoll = "future pregnancy intention";
run;
	
data mydata.LCA1119aim3;
set mydata.LCA1119aim3;

if best=1 then php=2;
else if best=2 then php=0;
else if best=3 then php=1;
if pregint=0 and certain =0 then intention=0; /* preg_intention 6 categorties*/
else if pregint=1 and certain=1 then intention=1;
else if pregint=1 and certain =0 then intention=2;
else if pregint=2 and certain=0 then intention=2;
else if pregint=0 and certain=1 then intention=3;
else if pregint=2 and certain=1 then intention=3;
else if pregint=1 and certain=2 then intention=4;
else if pregint=0 and certain=2 then intention=4;
else if pregint=2 and certain=2 then intention=5;

if intention=0 then intentcoll=0; /* preg_intention_four categories*/
else if intention=1 then intentcoll=1;
else if intention in (2:3) then intentcoll=2;
else if intention in (4:5) then intentcoll=3;
format php phpx. intentcoll intentcollx.;
run;
proc freq data=mydata.LCA1119aim3;
tables php intentcoll; run;
/*********************BIVARIATE ANALYSIS EXAMINING POSTIVE ASSOCIATIONS BETWEEN SDOH'S AND PREGINT, CERTAIN AND INTENTCOLL***************/

ods rtf file ="C:\Users\anakwea\Desktop\Men\RTF_results\2011_2019\Aim3\descriptives_by_intentcoll";
proc surveyfreq data=mydata.LCA1119aim3;
	tables  (php metro2 insur2 reprohlth workcons2 pir0 educ hisprace0 immi
	union2 UNION3 kiddos ages sexage2)*intentcoll / col row nowt cl chisq;
	tables (php metro2 insur2 reprohlth workcons2 pir0 educ hisprace0 immi
	union2)*pregint / col row nowt cl chisq;
	tables  (php metro2 insur2 reprohlth workcons2 pir0 educ hisprace0 immi
	union2)*certain /col row nowt cl chisq;
weight weightvar;
cluster panelvar;
strata stratvar;
run;
ods rtf close;

/********************FITTING A MULTINOMIAL REGRESSION MODEL FOR FERTILITY INTENTION AND PRECONCEPTION HEALTH*************/

ods rtf file="C:\Users\anakwea\Desktop\Men\RTF_results\2011_2019\Aim3\Crude Regressions";
proc surveylogistic data=mydata.LCA1119aim3; /*score test <.0001, therefore not an ordinal. Fit Multinomial regression*/
class php (ref="lowest risk") /param=ref;
model intentcoll(order=internal) = php; 
run;
/******** fitting the crude regression model "Definite NO" as ref********/
proc surveylogistic data=mydata.LCA1119aim3;
class php (ref="lowest risk")/ param=ref;
model intentcoll(ref="definite no") = php / link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;
/******** fitting the adjusted regression model "Definite NO" as ref - adjusted model********/

proc surveylogistic data=mydata.LCA1119aim3;
title "Fully adjusted model";
class php (ref="lowest risk") metro2 (ref='0')insur2 (ref="0") reprohlth (ref="2")
workcons2 (ref="2")pir0 (ref="1") educ0 (ref="1") hisprace2(ref="1")
immi(ref="0") union3(ref="0") kiddos(ref="0") ages(ref="1") sexage2(ref="1")/ param=ref;
model intentcoll(ref="definite no") = php metro2 insur2 reprohlth workcons2 pir0 educ0 hisprace0 immi
	UNION3 kiddos ages sexage2 / link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;

/******** fitting the adjusted regression model "Definite YES" as ref - CRUDE model********/

proc surveylogistic data=mydata.LCA1119aim3;
Title "Crude model, Definite yes as ref";
class php (ref="lowest risk")/ param=ref;
model intentcoll(ref="definite yes") = php / link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;

/******** fitting the adjusted regression model "Definite YES" as ref - Adjusted model********/
proc surveylogistic data=mydata.LCA1119aim3;
class php (ref="lowest risk") metro2 (ref='0')insur2 (ref="0") reprohlth (ref="2")
workcons2 (ref="2") workcons (ref="2") pir0 (ref="1") educ0 (ref="1") hisprace2(ref="1")
immi(ref="0") union3(ref="0") kiddos(ref="0") ages(ref="1") sexage2(ref="1")/ param=ref;
model intentcoll(ref="definite yes") = php metro2 insur2 reprohlth workcons2 pir0 educ0 hisprace0 immi
	UNION3 kiddos ages sexage2;  /*link=glogit; */
/*weight weightvar;
cluster panelvar;
strata stratvar;*/
run;


/******** fitting the adjusted regression model "Definite YES" as ref - drop metro, reprohlth, immigration, sexage********/

proc surveylogistic data=mydata.LCA1119aim3;
Title "Removing Metro, Reprohlth, Immigration, sexage from the model";
class php (ref="lowest risk") insur2 (ref="0") 
workcons2 (ref="2") workcons (ref="2") pir0 (ref="1") educ0 (ref="1") hisprace2(ref="1")
union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=ref;
model intentcoll(ref="definite yes") = php insur2  workcons2 pir0 educ0 hisprace0 
	UNION3 kiddos ages; /*link=glogit; */
/*weight weightvar;
cluster panelvar;
strata stratvar;*/
run;
/******** fitting the adjusted regression model "Definite NO" as ref - drop metro, reprohlth, immigration, sexage********/

proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting model with definite NO as reference";
class php (ref="lowest risk") insur2 (ref="0") 
workcons2 (ref="2") workcons (ref="2") pir0 (ref="1") educ0 (ref="1") hisprace2(ref="1")
union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=ref;
model intentcoll(ref="definite no") = php insur2  workcons2 pir0 educ0 hisprace0 
	UNION3 kiddos ages / link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;

/********************Examining the effcet of Far/contingent as baseline************************/
proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting model with definite NO as reference";
class php (ref="lowest risk") insur2 (ref="0") 
workcons2 (ref="2") workcons (ref="2") pir0 (ref="1") educ0 (ref="1") hisprace2(ref="1")
union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=ref;
model intentcoll(ref="Far or contingent intention") = php insur2 workcons2 pir0 educ0 hisprace0 UNION3 kiddos ages / link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;

ods rtf close;


/**************************** FITTING INTERACTION MODELS***********************************************/
/********MODEL 1: Interactive effect of RACE/ETHNICITY*********************************/
/********* Intentcoll= php*******************/

proc surveyfreq data=mydata.LCA1119aim3;
tables php intentcoll; run;
/*******************************************CRUDE intercation Model*************************************/
proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting model with definite NO as reference";
class  php (ref="lowest risk") hisprace2(ref="1")/ param=glm;
model intentcoll(ref="definite no") = php*hisprace2 / link=glogit expb; 
oddsratio php / at(hisprace2=all);
weight weightvar;
cluster panelvar;
strata stratvar;
run;

proc export data=mydata.LCA1119aim3 outfile="C:\Users\anakwea\Desktop\Men\AIM_3_stata\LCA1119aim3.dta"; run;

proc freq data=mydata.LCA1119aim3; tables hisprace0 hisprace2; run;
/****************************************************************************************************************************************
*****************************************************************************************************************************************
*****************************************************************************************************************************************
**********************************AIM 3B ANALYSIS****************************************************************************************/

/*********************************INTERACTION ANALYSIS (adjusted)**************************************************/
proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting intercation model for php*hisprace2";
class  php (ref="lowest risk") insur2 (ref="0") workcons2 (ref="2") pir0 (ref="1") educ0 (ref="1") 
hisprace2(ref="1") union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=glm;
model intentcoll(ref="definite yes") = php|hisprace2 insur2 workcons2 pir0 educ0 UNION3 kiddos ages/ link=glogit; 
lsmeans php*hisprace2 /oddsratio adj=bon cl diff;
slice php*hisprace2 /sliceby=hisprace2;
weight weightvar;
cluster panelvar;
strata stratvar;
run;

proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting intercation model for php*INSUR2";
class  php (ref="lowest risk") insur2 (ref="0") workcons2 (ref="2") pir0 (ref="1") educ0 (ref="1") 
hisprace2(ref="1") union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=glm;
model intentcoll(ref="definite yes") = php|insur2 hisprace2 workcons2 pir0 educ0 UNION3 kiddos ages/ link=glogit; 
lsmeans php*INSUR2 /oddsratio adj=bon cl diff;
slice php*INSUR2 /sliceby=INSUR2;
weight weightvar;
cluster panelvar;
strata stratvar;
run;

proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting intercation model for php*WORKCONS2";
class  php (ref="lowest risk") insur2 (ref="0") workcons2 (ref="2") pir0 (ref="1") educ0 (ref="1") 
hisprace2(ref="1") union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=glm;
model intentcoll(ref="definite yes") = php|workcons2 insur2 hisprace2  pir0 educ0 UNION3 kiddos ages/ link=glogit; 
lsmeans php*workcons2 /oddsratio adj=bon cl diff;
slice php*workcons2/sliceby=workcons2;
weight weightvar;
cluster panelvar;
strata stratvar;
run;

proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting intercation model for php*PIR";
class  php (ref="lowest risk") insur2 (ref="0") workcons2 (ref="2") pir0 (ref="1") educ0 (ref="1") 
hisprace2(ref="1") union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=glm;
model intentcoll(ref="definite yes") = php|pir0 workcons2 insur2 hisprace2 educ0 UNION3 kiddos ages/ link=glogit; 
lsmeans php*pir0/oddsratio adj=bon cl diff;
slice php*pir0/sliceby=pir0;
weight weightvar;
cluster panelvar;
strata stratvar;
run;

proc freq data=mydata.LCA1119aim3; tables educ0 educ education; run;

proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting intercation model for php*EDUC";
class  php (ref="lowest risk") insur2 (ref="0") workcons2 (ref="2") pir0 (ref="1") education (ref="0") 
hisprace2(ref="1") union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=glm;
model intentcoll(ref="definite yes") = php|education pir0 workcons2 insur2 hisprace2 UNION3 kiddos ages/link=glogit; 
lsmeans php*education/oddsratio adj=bon cl diff;
slice php*education/sliceby=education;
weight weightvar;
cluster panelvar;
strata stratvar;
run;

/*********************************STRATIFYING MODELS BY RACE: INTERACTION ANALYSIS (crude and adjusted)********************************************/

ods rtf file="C:\Users\anakwea\Desktop\Men\RTF_results\2011_2019\Aim3\Race_stratified";


proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting stratified model for RACE-crude";
class  php (ref="lowest risk")/ param=ref;
	DOMAIN hisprace2;
	model intentcoll(ref="definite no") = php /link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;

proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting stratified model for RACE-crude (substance use as ref)";
class  php (ref="lowest risk")/ param=ref;
	DOMAIN hisprace2;
	model intentcoll(ref="definite yes") = php /link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;

proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting stratified model for RACE";
class  php (ref="lowest risk") insur2 (ref="0") workcons2 (ref="2") pir0 (ref="1") educ (ref="0") 
hisprace2(ref="1") union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=ref;
	DOMAIN hisprace2;
	model intentcoll(ref="definite no") = php educ pir0 workcons2 insur2 UNION3 kiddos ages/link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;
proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting stratified model for RACE-adjusted(substance use as ref)";
class  php (ref="lowest risk") insur2 (ref="0") workcons2 (ref="2") pir0 (ref="1") educ (ref="0") 
hisprace2(ref="1") union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=ref;
	DOMAIN hisprace2;
	model intentcoll(ref="definite yes") = php educ pir0 workcons2 insur2 UNION3 kiddos ages /link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;
proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting stratified model for RACE-adjusted(substance use as ref)";
class  php (ref="lowest risk") insur2 (ref="0") workcons2 (ref="2") pir0 (ref="1") educ (ref="0") 
hisprace2(ref="1") union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=ref;
	DOMAIN hisprace2;
	model intentcoll(ref="Far or contingent intention") = php educ pir0 workcons2 insur2 UNION3 kiddos ages /link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;

ods rtf close;

ods rtf file="C:\Users\anakwea\Desktop\Men\RTF_results\2011_2019\Aim3\Educ_stratified";
/*********************************STRATIFYING MODELS BY EDUC STATUS: INTERACTION ANALYSIS (adjusted)*******************************/
/*****creating a new " education variable with 3 categories************/
data mydata.LCA1119aim3;
set mydata.LCA1119aim3;
if educ in (0:1) then education=0; else if educ=2 then education=1; else if educ=3 then education=2;
run;
proc freq data=mydata.LCA1119aim3; tables educ education educ*education; run;


proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting stratified model for Educ";
class  php (ref="lowest risk")/ param=ref;
	DOMAIN education;
	model intentcoll(ref="definite no") = php pir0 workcons2 insur2 UNION3 kiddos ages/link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;
proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting stratified model for Educ -crude (substance use as ref)";
class  php (ref="lowest risk")/ param=ref;
	DOMAIN education;
	model intentcoll(ref="definite yes") = php pir0 workcons2 insur2 UNION3 kiddos ages/link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;

proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting stratified model for Educ";
class  php (ref="lowest risk") insur2 (ref="0") workcons2 (ref="2") pir0 (ref="1")  
hisprace2(ref="1") union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=ref;
	DOMAIN education;
	model intentcoll(ref="definite no") = php pir0 workcons2 insur2 UNION3 kiddos ages/link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;
proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting stratified model for Educ adjusted (substance use as ref)";
class  php (ref="lowest risk") insur2 (ref="0") workcons2 (ref="2") pir0 (ref="1") 
hisprace2(ref="1") union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=ref;
	DOMAIN education;
	model intentcoll(ref="definite yes") = php pir0 workcons2 insur2 UNION3 kiddos ages/link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;
ods rtf close;

/**********************Forest plots for Race stratified models******************************************/
ods output "odds ratios" = raceforest;

proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting stratified model for RACE";
class  php (ref="lowest risk") insur2 (ref="0") workcons2 (ref="2") pir0 (ref="1") educ (ref="0") 
hisprace2(ref="1") union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=ref;
	DOMAIN hisprace2;
	model intentcoll(ref="definite yes") = php educ pir0 workcons2 insur2 UNION3 kiddos ages/link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;
/*proc surveylogistic data=mydata.LCA1119aim3;
Title "Fitting stratified model for RACE-adjusted(substance use as ref)";
class  php (ref="lowest risk") insur2 (ref="0") workcons2 (ref="2") pir0 (ref="1") educ (ref="0") 
hisprace2(ref="1") union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=ref;
	DOMAIN hisprace2;
	model intentcoll(ref="definite yes") = php php educ pir0 workcons2 insur2 UNION3 kiddos ages /link=glogit; 
weight weightvar;
cluster panelvar;
strata stratvar;
run;

/**************proc print forest plots*********************************************************/
data raceforest; set raceforest;
effect=upcase(effect); run;
Title "Forest plot for Race stratified models";
proc sgplot data=raceforest;
	scatter x=oddsratioest y=effect/ xerrorlower=lowercl xerrorupper=uppercl markerattrs=or
	(symbol=DiamondFilled size=8);
		refline 1/ axis=x;
	xaxis label ='OR and 95% CI' min=0.0 max=2.5;
	yaxis label ='Covariates';
run;


proc freq data=mydata.lca1119aim3; tables age ages; run;

proc freq data=lca1119aim3; tables P12MOCON; run; 





/**********************MULTINOMIAL POWER ANALYSIS********************************/
proc logistic data=mydata.LCA1119aim3;
class php (ref="lowest risk") metro2 (ref='0')insur2 (ref="0") reprohlth (ref="2")
workcons2 (ref="2") workcons (ref="2") pir0 (ref="1") educ0 (ref="1") hisprace2(ref="1")
immi(ref="0") union3(ref="0") kiddos(ref="0") ages(ref="1") sexage2(ref="1")/ param=ref;
model intentcoll(ref="definite yes") = php metro2 insur2 reprohlth workcons2 pir0 educ0 hisprace0 immi
	UNION3 kiddos ages sexage2/ link=glogit; 
run;


/******** fitting the adjusted regression model "Definite YES" as ref - drop metro, reprohlth, immigration, sexage********/

proc logistic data=mydata.LCA1119aim3;
Title "Removing Metro, Reprohlth, Immigration, sexage from the model";
class php (ref="lowest risk") insur2 (ref="0") 
workcons2 (ref="2") workcons (ref="2") pir0 (ref="1") educ0 (ref="1") hisprace2(ref="1")
union3(ref="0") kiddos(ref="0") ages(ref="1")/ param=ref;
model intentcoll(ref="definite yes") = php insur2  workcons2 pir0 educ0 hisprace0 
	UNION3 kiddos ages/ link=glogit; 
run;
/******POWER ANALYSIS FOR MULTINOMIAL REGRESSON************************/
proc power;
custom
dist = chisquare 
primnc = 0.00373765 /*-2LL(full)- (-2LL (reduced))/df(reduced)*/
testdf = 33
ntotal = .
alpha = 0.05
power = 0.90;
run;


PROC SGPANEL DATA=mydata.LCA1119aim3 NOAUTOLEGEND;
	PANELBY hisprace2 / ROWS=1 COLUMNS=5 NOVARNAME SORT=DATA NOHEADERBORDER HEADERBACKCOLOR=WHITE HEADERATTRS=(FAMILY = ARIAL SIZE = 12 WEIGHT = BOLD);
	WHERE OUTCOME1 = 'intentcoll';
	REFLINE 0 / AXIS=X LINEATTRS=(THICKNESS = 3 COLOR = DARKGRAY PATTERN = DASH);
	SCATTER X=ESTIMATE Y=effect / XERRORLOWER=LCL XERRORUPPER=UCL NOERRORCAPS ERRORBARATTRS=(THICKNESS=3) MARKERATTRS=(SYMBOL = CIRCLEFILLED SIZE = 15)
		COLORMODEL=(BLACK PURPLE DARKRED DARKGREEN DARKBLUE) COLORRESPONSE=ORDER;
	COLAXIS LABEL='Effect on intentcoll' LABELATTRS=(FAMILY = ARIAL SIZE = 12 WEIGHT = BOLD) VALUEATTRS=(FAMILY = ARIAL SIZE = 10 WEIGHT = BOLD);
	ROWAXIS LABEL='effect' LABELATTRS=(FAMILY = ARIAL SIZE = 12 WEIGHT = BOLD) VALUEATTRS=(FAMILY = ARIAL SIZE = 10 WEIGHT = BOLD);
	TITLE 'Effect of PHP on fertiity intentions';
RUN;

proc freq data =mydata.LCA1119aim3;
tables hisprace2; run;
